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Abstract 

Numerical simulations are used to test the kinetic theory constitutive relations of inertial granular 
shear flow. These predictions are shown to be accurate in the dilute regime, where only binary 
collisions are relevant, but underestimate the measured value in the dense regime, where force 
networks of size £ are present. The discrepancy in the dense regime is due to non-collisional 
forces that we measure directly in our simulations and arise from elastic deformations of the force 
networks. We model the non-collisional stress by summing over all paths that elastic waves travel 
through force networks. This results in an analytical theory that successfully predicts the stress 
tensor over the entire inertial regime without any adjustable parameters. 
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I. INTRODUCTION 



Granular materials exhibit a broad spectrum of behaviors that have been difficult to 
capture theoretically, especially in the dense regime. A fundamental question is when, and 
whether, collective motion becomes important for understanding the macroscopic state of 
the system l|. This is a particularly complex issue for granular flows, where the material 
structure remains amorphous on all scales. While a noticeable structural change would surely 
help pinpoint degrees of freedom that govern collective dynamics, granular flows require us 
to use the dynamics to search for structure. 

A central tool in this process is numerical simulation, and new insights from simulations 
can provide important clues for theory. In the companion paper j^] we have presented 
evidence for a growing length scale related to correlations between grain forces. These 
correlations arise from force chain networks that span the space between grains and transmit 
forces elastically. The emergence of force networks alters the mechanisms of momentum 
transfer, and contact force statistics depend on the size £ of force networks. 

Since contact forces are sensitive to the value of £, the derivation of constitutive relations 
to describe the stress tensor must also be based on properties of the networks. While network 
properties have been incorporated into previous studies of static granular packings 
a, |7|, |8|, |9|, the importance of force networks in the inertial regime has not been adequately 
taken into account. Instead, constitutive relations are generally obtained using comparisons 



with liquids [10 



dense materials 



111 



and especially kinetic theory 12 



16 



which can be extended to treat 



17l | as long as correlations beyond nearest neighbors are absent. 



In this paper we examine two models for constitutive relations in granular flow: kinetic 
theory and a new model, referred to as the force network model. Kinetic theory is rooted in 
the dilute limit and assumes no long range correlations between grains. The force network 
model is motivated by our simulations and attempts to capture the effects of finite sized 
force networks in a mean-field framework, which spans system behavior from dilute to dense 
regimes. 

The predictions of these models are compared using measurements of the stress tensor 
E Q( 3, which is given by 

N c 



E Q/3 y = ^m*K- MQ )(^- M/3 )+ 

i=l {i,j}=l 



(1) 
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In this equation, Greek subscripts denote components, Roman superscripts denote grains, 
and V represents the volume. The first term quantifies stress resulting from fluctuating 
velocities, which is related to the mass of each grain m % and the difference between grain 
velocity v J and average velocity u (where bold face symbols denote full, three component 
vectors). The second term arises from contacts between grains: it depends on the contact 
forces F lJ and the distance between contacting grains a %3 . The sum is taken over all contacts 
{i, j} in the system. 

We focus on the second term, called the static stress, since it is directly associated with 
contact forces, and we specialize to non-frictional systems where cr* J x F 13 = 0. This allows 
us to write the static stress tensor E^g as 

c 

K,V= Yl (2) 
{m}=i 

where a 13 = a %3 a %3 and F %3 is the magnitude of the contact force between pairs {i, j} of 
contacting grains. 

We begin in Section[TTlwith an investigation of dilute inertial flows and test the predictions 
of hard-sphere kinetic theory. In Section II III we focus on dense inertial flows and introduce 
the force network model. 



II. DILUTE INERTIAL FLOWS 

Over the past twenty five years 18 



191 ] the kinetic theory of dense gases has been general- 
ized to include granular flows, where thermal fluctuations are absent and energy is dissipated 
at each contact. The dissipation mechanism most often considered is instantaneous collisions 



with constant restitution coefficient- this is called hard-sphere kinetic theory 2Q( ■ In order 
to make progress using hard-sphere kinetic theory, it is necessary to begin by postulating that 
only binary collisions occur between grains. Without this assumption, calculations quickly 
become intractable since high order correlations must be included in kinetic integrals. 

The principal microscopic input to hard-sphere kinetic theory is the collision rule between 
grains. This relates the initial velocities of two interacting grains {v n , V 3 } to their final 
velocities {v l , v 3 }: 

( v i _ v i) . &io = _ e ( v " _ v '*) . a i3 . (3) 
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The normal coefficient of restitution e that appears in this equation determines the energy 
dissipated in each collision: for e = 1 the system is elastic and no energy is dissipated; as e 
is reduced to zero the energy dissipation scales as 1 — e 2 . 

Kinetic theory relies on the assumption that only binary collision are relevant and is 



therefore expected to 



tacts arise 
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jreak down as the density of the flow increases and long-lasting con- 



24j . Quantitative bounds over which the binary collision assumption 



holds have only recently been estimated 
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26l |. The Contact Dynamics (CD) algorithm 



used here to simulate granular flows is well suited for testing the relevance of the binary 
collision assumption and bounding the dilute regime. Like hard-sphere kinetic theory, the 
CD algorithm employs a normal coefficient of restitution. However the CD algorithm does 
not assume a priori that only binary collisions occur. 

On the contrary, we observe that multi-grain contacts are the dominant interaction in 

n 

dense flows. This is evident in force correlation measurements [2j, which identify a growing 
correlation length £. This length scale should be an indicator of the breakdown of the binary 
collision assumption since it implies that grain forces are correlated over large distances and 
do not simply depend on nearest neighbors. To see how this comes about, it is useful to 
measure the static stress tensor. 

Equation ([2]) gives the microscopic expression for the static stress tensor, which depends 
on contact forces between grains. In the case that only binary collisions are considered, as 
in kinetic theory, the value of the static stress tensor is determined by inserting the binary 
collisional force into Equation (j2J). The collisional force is the force that occurs between a 
pair of colliding grains that are not part of a force network. Given the initial velocities V 1 
of grains i, the collisional force between two grains is given by 

Fli = (1 + e)^ [(v' J - v") ■ *«] /dt, (4) 

where e is the normal restitution coefficient, /i = m % m> j{m % + m?) is the reduced mass, and 
& tJ is the unit vector connecting the centers of grains i and j. This expression is equal 
to the change of momentum of grains % and j, per simulation time step dt, due to binary 
interactions. Since all of the parameters can be measured in simulations, the collisional 
force is a useful probe of the dynamics and was used in the companion paper to exhibit the 
presence of force networks. 
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Inserting the collisional force into Equation (j2j) yields the "collisional" stress tensor 

c 

= £ ^> iiF t (5) 

0,i}=o 



which is the stress tensor that all hard-sphere kinetic theories attempt to model [26|]. This 
is because only binary collisions between grains are considered and contact forces can be no 
larger than the value given in Equation (pE|). 

If relevant interactions occur only via binary collisions, the collisional stress in Equa- 
tion ([5]) and the static stress in Equation ([2]) are identical. However, a discrepancy between 
the collisional and static stress tensors indicates that momentum transfer is not carried 
out by binary collisions alone. Comparing the static and collisional stress tensors in Equa- 
tions (^j and §5§ therefore provides an opportunity to test the binary collision assumption 
and thereby determine when kinetic theory can be applied to hard-sphere granular flows. 

In Figures Q] and [2] we plot measurements of the static and collisional stress tensor, for a 
wide range of restitution coefficients and packing fractions, in terms of the pressure p and 
shear stress s, which are related to the stress tensors by 

{P s ,p hc } = {^(S s n + £ s 22 ), + ^ 2 C )}, (6) 

{s s ,s hc } = {S S 12 ,S^} = {S & 21 ,S^}. (7) 

In Figures Q] and [2] both the collisional and static values of the pressure and shear stress 
are normalized by common factors that are explained later, in the paragraph beneath Equa- 
tion (Tbfl) . For now, it is only important to note that there is a regime where the binary 
collision assumption holds and the normalized stress tensors are equal. The bounds of this 
regime depend on the value of both the restitution coefficient and the packing fraction. The 
data supports a conclusion that the dilute regime is approached as density is reduced or 
restitution is increased. 

Instead of characterizing the dilute regime in terms of restitution and packing, it is useful 
to relate it to the length scale £. In Figures [1] and [2] we have colored the data points where 
£/£ei > 1.25. For all of our data, this simple condition on £ nicely characterizes the dilute 
regime- if £/£ e i < 1-25 then the static stress tensor is approximately equal to the collisional 
stress tensor and the predictions of kinetic theory apply; if £/£ e i > 1-25 then interactions 
between networks of grains become important and kinetic theory does not apply. The 
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crossover value of £/£ e i = 1-25 is also where the contact force distribution P(f) loses its 
peak [2]. This provides further evidence that £/£ e i = 1-25 is a good quantitative bound for 
the dilute regime. 

In the dilute regime where £/£ e i < 1-25 and Eig = S^, we can test the predictions 
of hard-sphere kinetic theory for _both__the pressure and shear stress. These predictions 



have been obtained recently lis, using the Chapman- Enskog expansion to solve 
the Enskog equation. The Enskog equation determines the time dependence of the one- 
particle probability distribution function (pdf) in terms of collision events between grains. 
Collision events consist of binary interactions and the time dependence of the one-particle 
pdf can therefore be expressed in terms of just the two-particle pdf. For hard-sphere granular 
materials, Enskog's equation reads 

(dt + vx ■ VxJ/^Crx.Vx.t) = J^Vx], (8) 

with Je given by 

■Mrx,vx] = a J dv 2 J d&B(& ■ g)(<x • g) (9) 
x[e- 2 / (2) (ri,r!- <r,v'x,v^) 
-/ (2) (r 1 ,r 1 + < T,v 1 ,v 2 ,t)]. 

Physically, the time dependence of the one-particle pdf /W is related to a collisional term Je 
that quantifies the probability to gain and lose contributions at a certain velocity Vx- The 
first term in Je gives the probability that a binary collision between grains results in a grain 
having velocity v x , and the second term gives the probability that a binary collision occurs 
involving a grain that has velocity v 1; thereby reducing /^(vx). The binary collisions occur 
according to Equation ([3]) and primed velocities represent pre-collisional values. G is the 
step function and g = v x — v 2 . For hard sphere granular flows, this form for the Enskog 
equation can be formally derived, starting with the binary collision assumption and the 



pseudo-Liouville equation [2CJ, |27|] . 



A prediction for the collisional stress tensor is obtained by multiplying Equation (jSJ) on 
each side by mv 1; where m is the particle mass, and integrating over v x . This yields the 
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FIG. 1: (Color online). Main Figures: Normalized values of the pressure (left) and shear stress 
(right) for large restitution coefficients e = 0.92, e = 0.75 and e = 0.5. The pressure is normalized 
by po from Equation (|12p and the shear stress is normalized by jrjo from Equation (|13p . The dilute 
regime is characterized by the range of restitution and packing where the static and collisional values 
are equal. Filled data points correspond to values of restitution and packing where £/£ei > 1.25- 
this provides a simple quantitative condition for the boundary of the dilute regime. Kinetic theory 
is expected to apply in the dilute regime and the prediction for the pressure is accurate for all 
e. The prediction for the shear stress overestimates the actual value, due to positive velocity 
correlations. Insets: Pre-collisional velocity correlations as a function of packing fraction (defined 
in Equation (fT5"|) ). 




FIG. 2: (Color online). The same as Figure (pQ), for small restitution coefficients of e = 0.25, 
e = 0.1 and e = 0. 



transport equation for momentum density 151 ]. which gives the stress tensor 



1 + e 



ma J dvi J d\2 J d&Q(& ■ g)(<x ■ g,) 2 cr a crp 



/ d\fW [r - (1 - A)cr, r + Act, v 1; v 2 , t] . 
Jo 



(10) 



In order to determine the stress tensor and solve the Enskog equation, it is necessary to 
express in terms of Assuming there are no velocity correlations between grains 



S 



that are about to collide yields 

f [2] (ri, r 2 , Vl , v 2 , t) = X (ri, r 2 )/(r ls v ls t)f(r 2 , v 2 , t) (11) 

and reduces the Enskog Equation (jSJ) to a non-linear differential equation for the one-particle 
pdf. The function \ is interpreted as the equilibrium correlation function at contact and 
depends on the local value of the density. 

Once the Enskog equation has been expressed in terms of only the one-particle pdf, it 



can be solved using the Chapman-Enskog expansion 28], |29J, which expands f^ 1 ' and Je 



in gradients of the mass density, momentum density, and energy density. This pro cess has 



been carried out for granular shear flows to first order in the gradients in Refs 



In two dimensions, this gives predictions 13] for the pressure pP red and shear stress sP red 



. 3, HQ- 



„pred 

= (1 + e) X u, (12) 



Po 

s prcd 4^/4 



+ vx (l + e)f{e)), (13) 



7770 vr \5 — e 

,l ) 5-e \ 81 - lie + :M)eHl - e) ) ' ( ' 

which only depend on the restitution e, the packing fraction u, and the pair correlation func- 
tion at contact x- The normalizing factors are given by po = nmT/2 and rjo = m/a^T/2-K, 
where m is the grain mass, a the grain diameter, T the granular temperature, and n the 
number density. We measure all the variables of Equations (|12p and ([TBI in simulations and 
test the predictions of hard-sphere kinetic theory without fitting parameters. We use the 
average value of grain mass and diameter for m and a, and determine x by tracking the 
number of collisions that occur per second (which we denote by u) in equilibrium simulations 



where e = 1 and 7 = 0. Enskog theory relates x to oj through the equation uj = V2ii5Txna. 



30] 



This method for measuring x has been used in other recent studies 

We plot the normalized predictions from Equations (11211141) in Figures [1] and [21 where 
they are compared to data for the stress tensor. Since hard-sphere kinetic theory assumes 
binary collisions, the predictions only apply to the dilute regime where £/£ e i < 1-25. These 
are the open symbols in Figures [TJ and [2J We immediately notice that the prediction for the 
pressure matches the measured pressure in all of the dilute systems we have investigated. 
Considering that there are no adjustable parameters, this is a tremendous success for kinetic 
theory. 



The prediction for the shear stress matches the data in the dilute regime only for large 
restitution. As the restitution becomes smaller, the prediction for the shear stress begins 
to overestimate its measured value. This overestimation is due to pre-collisional velocity 



correlations, a mechanism that has been investigated in previous studies [31|- Since Equa- 
tion (1111) assumes no pre-collisional velocity correlations, if these correlations exist then the 
average momentum transferred in each collision changes. If the pre-collisional normal veloc- 
ities of two grains tend to be aligned (anti-aligned) then the average momentum transferred 
decreases (increases), causing the kinetic theory prediction to overestimate (underestimate) 
the measured values. 

In the insets of Figures [T] and [2J we plot measurements of the pre-collisional velocity 
correlations C v defined as 

C v = ((v n -cr^)(V> -&v))/T, (15) 

which are normalized by the granular temperature. This definition yields a positive value 
when pre-collisional grain velocities tend to be aligned, and for all restitution coefficients 
we observe that the correlations are positive. In addition, the magnitude of the discrepancy 
between measured and predicted shear stress is roughly proportional to the size of the 
velocity correlations. These observations support the conclusion that the errors in the kinetic 
theory prediction are due to pre-collisional velocity correlations. 

To summarize, the stress tensor predicted by hard-sphere kinetic theory matches data 
from our simulations in the dilute regime where £/£ e i < 1-25 and at high restitution coeffi- 
cients. As the restitution coefficient is decreased in dilute flows, the prediction for the shear 
stress begins to fail due to the "molecular chaos" assumption of Equation (fTTj) . Additionally, 
as £/£ei becomes larger than 1.25 in the dense regime, the hard-sphere kinetic theory predic- 
tions are inaccurate since the binary collision assumption is not valid. These observations 
indicate that the fundamental assumptions of molecular chaos and binary interaction do not 



always apply, and must 



Recent research 



20 



3e addressed. 



321 ] has concentrated on refining the molecular chaos assumption 



to account for velocity correlations in the dilute regime, which have been measured exten- 



sively 3l|, |33|, |3J, |35j, |36[]. However in the dense regime, where networks of interacting 



grains become important, even an exact inclusion of velocity correlations will not accurately 
describe the physics since the binary collision assumption does not apply. 
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III. DENSE INERTIAL FLOWS 



Dense inertial flows are not quasi-static and cannot be modeled by assuming binary 
interactions between grains. They span the range of densities between the dilute and quasi- 
static regimes, and exhibit properties of both limits. Like dilute flows, dense inertial flows 
are characterized by a Bagnold rheology where the stress tensor is proportional to the square 



of the shear rate 



37 



38 



39|. However, as in quasi-static flows, the value of the stress tensor 



also depends on the properties of force chain networks. In fact, the size of the networks is 
a parameter that is relevant in all regimes of granular flow. It varies continuously from the 
dilute regime, where interactions are binary and the network size is unity, to the quasi-static 
regime where force chain networks extend over the entire system. Since the behavior of 
granular flows can be nicely characterized in terms of the network size, we focus here on 
the role of force chain networks in determining the value of the stress tensor and construct 
a force network model of momentum transport in granular materials. The basic idea is 
that forces are transferred elastically through finite sized contact networks and the value of 
the contact force between any pair of grains depends on both the relative velocity of the 
contacting pair (a "collisional" contribution) and the values of the other contact forces in 
the network, even those a long distance away (an "elastic" contribution). This leads to a 
prediction for the stress tensor that holds over the entire inertial regime and incorporates 
the non-collisional stress that arises when force networks have formed. 



A. The Force Network Model 

A central feature of dense granular materials is that forces can be transferred over dis- 
tances much larger than the grain size. This is especially evident in static packings of grains, 
but is also important when considering dense inertial flows. The only necessary requirement 
for force propagation is the existence of connected networks of interacting grains. When 
this requirement is met, elastic waves propagate through the networks at a speed set by the 
values of the moduli and there is an elastic component of the stress tensor. This elastic 
contribution must be added to the collisional part of the stress in Equation (jSJ). Therefore 
the full static stress tensor is comprised of two terms: one describing the elastic response 
and one describing the collisional response. 
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From a microscopic viewpoint, the value of the stress tensor is determined by contact 
forces between grains, as in Equation ([2]). Therefore, in the presence of force networks, 
these contact forces must also be comprised of collisional and elastic terms. The collisional 
term is given by Equation (J3|) and represents the local force due to collisions between pairs 
of grains. It depends only on properties of the two contacting grains. The elastic term is a 
result of elastic deformations in the contact network. It is a non-local force that arises from 
the network applying an effective pressure on every pair of contacting grains. 

In dilute flows only the collisional part of contact forces is non-zero and the stress tensor 
is well described by kinetic theory. In the quasi-static regime the elastic part of the forces 
is much larger than the collisional part, and the latter is usually disregarded. In dense 
inertial flows, both terms contribute. The force network model extends hard-sphere kinetic 
theory by explicitly calculating the effects of elastic waves in force networks. This leads to 
a prediction for the elastic portion of contact forces, which result from forces propagating 
through force networks at the elastic wave speed. Inertial flows correspond to the limit where 
forces propagate throughout the entire network before it is destroyed. This is equivalent to 
the limit where the elastic wave speed is infinite, or the grains are perfectly rigid. 

Mathematically, the contact force F v between grains i and j is equal to the sum of a 
collisional term, plus elastic effects from the network: 



In this equation the first term is the collisional force, defined in Equation (J4j), and the 
second term arises from forces that propagate through paths in the force chain network, as 
illustrated in Figure El We split this term into contributions that represent added forces 
from different path-lengths £. The total elastic force created by the network is therefore 
equal to the sum of the contributions Jy 7 over all possible path-lengths £ < £ max in the force 
network. 

Figure [3] illustrates force propagation and defines the notion of path- length. For example, 
because grain d is in contact with grain e, this increases the force between grains e and /, 
which increases the force between grains / and i, which has the net effect of increasing the 
contact force between grains i and j. We denote this as a path of length three (£ = 3) since 
the local force must propagate through three links to influence the contact between % 
and j. 



4 




(16) 
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FIG. 3: A network of contacting grains. The contact force between the shaded grains i and j 
is determined by the local collisional force F?~ , plus elastic contributions from forces propagating 
through the network. It is convenient to organize these non-local effects into contributions from 
different path-lengths. On the top left is a path of size three where the local contact force between 
grains d and e is transferred through the network to grain i. On the top right is a path of size two, 
and on the bottom right is a path of size one. 

The net effect of forces propagating through different path-lengths can be calculated 
explicitly in the limit of large elastic wave speeds. We begin by considering a path of length 
I = 1, as illustrated in Figure [3j because grain a is in contact with grain j, the local contact 
force F^. increases the value of i™ by an amount equal to multiplied by the cosine of 
the angle between the unit vectors connecting contacts {a,j} and {i,j}- If we assume that 
grain i has Zj contacts labeled by m and grain j has Zj contacts labeled by n, then the effect 
of all paths of length one is to increase i™ in Equation f|T6|) by an amount 



where & a is the unit vector connecting the center of grains a and b. This expression includes 
all of the effects from paths of length I = 1 on each of the contacting grains i and j. 

In an analogous manner, the path from grain b to c to j in Figure [3] comprises a path of 
length two {£ = 2), which increases the value of F iJ due to the local force between b and c. 
Note that we ignore the fact that the local force between c and j also increases F % \ since 
this contribution was already included in Equation ffTTj) . The total additional force between 




(17) 
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grains i and j arising from paths of length two is given by 



Zm 

rrdj \ r ~ mi ~ ij \ r ~ pm ~ mi T^V m 



+ °- nj -° l ' J & qn -v nj F^, (18) 

where once again grain % has Zi contacts labeled by m and grain j has Zj contacts labeled by 
n. To calculate the effect of paths of length two, we also take into account the z m contacts 
of grain m, labeled by p, and the z n contacts of grain n, labeled by q. 

The contributions to from arbitrary path-lengths I can be determined by continuing 
the above arguments. They depend on the coordination number z and are also sensitive to 
the geometric arrangement of force networks. One important constraint arises as z becomes 
large, which is a result of energy conservation. If we consider the total energy T* J that an 
arbitrary pair of contacting grains % and j contribute to the network, this must always be 
less than or equal to the total kinetic energy initially stored in the contact (via the kinetic 
energies of grains i and j). This is because, while energy is conserved in this process, some 
is transferred to the elastic energy of the network and some remains as kinetic energy of 
grains i a nd j . The average collisional force (F^ J C ) is proportional to the square of the relative 



velocities 15( and thereby proportional to the kinetic energy of the contact. Therefore, on 
average, the magnitude of the collisional force transferred to nearest neighbors must not be 
greater than (F^ J C ). This leads to the constraint equation 

< (19) 

which is a global constraint. It ensures that the total energy supplied to the network never 
exceeds the initial kinetic energy as the elastic waves move from first nearest neighbors, to 
second nearest neighbors, and so on. 

Equations (11711191) model the physical origin of elastic forces that exist in dense granular 
materials and allow for a complete determination of J-f . The numerical value of the elastic 
force between a pair of contacting grains is calculated by summing these contributions over 
all possible path- lengths £, as in Equation ffTBT) . 

The maximum possible path length £ max is constrained by the size of the force networks. 
A straight chain of grains that spans the network has £* — £/£ e i — 1- While there are also 
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network spanning chains with £ > £*, their contributions to elastic forces are diminished 
since the magnitude of the collisional force at the end of the chain is reduced by the product 
of the cosine of the angles in the chain. We will therefore set £ max = £/£ei — 1 and only 
consider the network spanning chains that give the largest contribution. This amounts to 
completely ignoring the effects of path lengths with I > £ max which increase the elastic force 
on the contact {i, j} by an amount A = ^^ max+1 J~\ 3 ■ We do not expect this approximation 
to produce a large error in the total elastic force on {i,j} since A < and, in general, 

<F £+1 ^ F \ for £ < £ max . 

Given the above analysis, it is possible to completely determine the stress tensor based 
on properties of the force network and the collisional stresses. This is carried out in Subsec- 
tion HIIE] and the predictions are tested in Subsection IIII CI 

B. Calculating the stress tensor 

We calculate the stress tensor by rewriting the equations of the force network model in 
terms of integrals instead of sums. Here we carry out this substitution for a two-dimensional 
system, although it can be generalized to higher dimensions. If we consider the average force 
between two grains that contact at an angle 9, denoted F(9), then Equation (fl6l) can be 
generalized to read 

F{9) = F bc (9) + Ftf)- ( 2 °) 
t=\ 

This equates the average force between grains contacting at angle 9 to the average collisional 
force at that angle, plus elastic effects from the network. In what follows we will measure 
9 with respect to the axis of shear, so that 9 = corresponds to the axis where there is no 
gradient in grain velocities. 

Generalizing Equations (|17[) and (|18|) to arbitrary path-length £ gives 

l / .6» i _ 1 +27r/3 

F e (9 ) = II / d9 i cos(9 t -9 l _ 1 )C(9 t )F bc (0e) 

x (2i) 

In this equation, the sums from the Equations (TTTT) and (fI8]) have been replaced by integrals 
over 9i, which is the angular orientation of each link in the chain. Each integral contains 
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a cosine that replaces the dot product. Note that the bounds of the integrals are arranged 
such that the grain at link i is not permitted to overlap the grain at link i — The function 
Fbc{9e) provides the collisional force at the end of the path. In order to properly characterize 



the probability to have a contact at 0$, we also introduce the functions z 



C{9) gives the probability to have a single contact at angle 9 40, 



41 
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C 



43 



and P. 



Z (ey 



44J. In the 



case that there are two (or more) contacts on a single grain, which is necessary to form a 



chain, this probability must be modified 



45] . The function C(9 i )P z{ e i) (9 i - f _ a ) gives the 



conditional probability to have a contact at 0j, provided that there already is a contact at 
0j_i and there are a total of z{9i) contacts on the grain at 9i. The function p(z(9i)) is the 
probability to have a grain with z{9j) contacts, and the term in the parentheses captures 
the average effect of all possible combinations of contact numbers and angles on each grain 
in the force chain. 

It is useful to simplify Equation (I2T!) by defining the function Pc{9i — 0i-i), where 
C(9i)Pc{9i — 0j-i) gives the probability to have a contact at 0j, given that a contact al- 
ready exists at 0j_i, averaged over all possible fluctuations in z{9i). If we define 

(z - l)P c {9i - 6i- x ) = $>(^)) " l)iW0< - Oi-i), (22) 



t /•6» i _ 1 +27r/3 

Fi(0o) = T[(z ~ 1) / Mi cos(0, - 9 i „ 1 )P c (9 i - 9 l _ l )C(9 l )F hc (0i). (23) 



where z is the average coordination number in the entire network, then Equation (1211 
simplifies to 

f9i_i+27r/3 

Equation fl2"3"|) generalizes the sums in Equations f[T7|) and ffTSl) to arbitrary path length I. 
It does so by averaging over one-particle distribution functions. These come in the form 
of the average collisional force Fb c {9), the average contact probability on a single grain 
C(9i)Pc{9i — 0i_i), and the average coordination number z. Naturally, there must not be 
correlations between these variables in order for the predictions to apply. The definition of Pc 
ensures that there are no correlations between the contact probabilities and the coordination 
number, but it is not guaranteed that the collisional forces are not correlated with the 
network structure. In fact, when closed loops form in the force networks, the collisional 
forces and network structure do become correlated. Therefore, Equation ff23l only applies 
as a mean-field approximation that ignores the correlations induced from loops in the force 
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networks. This is a valid approximation since the effect of a collisional force that propagates 
around a loop is reduced by the product of the cosine of the angles around the loop, which 
is quite small compared with direct propagation. 

To solve the integrals in Equation (|23p . we change integration variables to X{ = 9,i — 6i-\. 
This results in the expression 



v. 



i=l 



2n/3 



2vr/3 



dxi cos(xi)P(xi)C 



3=1 



<>+£ 



(24) 



which incorporates the average effect of all forces that propagate through paths of length t 
on contacts with orientation 9. 

In addition to force propagation based on the cosine of the angle between subsequent 
contacts, we must also incorporate the global constraint from Equation (fi~9l) . This can be 
generalized to 



2tt/3 

T(6) = (F hc (0))(z-l) I dxP c (x)cos(x)C(6 + x)<(F hc (0)), 

-2vr/3 



(25) 



and restricts the total energy transferred through the network. 

Equations fl24|) and fl25|) . combined with the basic force network Equation fl20l) . comprise 
the integral form of the force network model. In order to carry out the integrations, it is 
necessary to know the functional form of C{9) and F^dO). These functions are ^-periodic 
and can be written as Fourier Series, keeping only terms that are also 7r-periodic. Previous 



research on the contact probability 
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41 



42 



43 



44j has shown that C{9) is well approxi- 



mated by keeping only the lowest Fourier terms. We find that F^ c {9) has the same property. 
We therefore approximate 

1 



C{9) 



2tt 



'1 + a c sin 26 + a' c cos 29), 



F hc (6) = (F hc )(l + a f sm29 + a' f cos26). 



(26) 
(27) 



In Figure H] we plot data of these functions for a granular material with e = and v = 0.79, 
along with a fit to the above equations. The fit is constructed by computing the fabric 
tensor <fi a p = {a a ap) and force-fabric tensor Xa/3 = {FhcVa<tp) / (Fbc) , where a a is the a 
component of the unit vector between contacting grain centers, F^ c is the collisional force 
on the contact, and the average is taken over all contacts. The anisotropies in the contact 
and force distributions {a c ,a' c ,af,a'f} are related to eigenvalues of the fabric tensors [431 ]. 
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FIG. 4: Polar plots of measurements (data points) and fits (lines) of (a) the contact probability 
distribution C{9) and (b) the collisional force distribution -Fbc(^)/(^bc) for a granular material 
with e = and v = 0.79. The lines are fits of Equations (|26p and (|27p . and the values of the 
Fourier components are plotted in Figure [5] 

which are simple to measure. We see from the plots that this first-order approximation for 
the contact probability and collisional force is quite good. 

We have measured C{6) and F^ c {6) for a wide range of restitution coefficients and packing 
fractions. In Figure [5] we plot the value of the Fourier components from Equations ([26]) 
and (]2Tj) . which characterize the functional form in all cases. These plots reveal that the 
anisotropy in both the contact probability and collisional force depends sensitively on the 
value of the packing fraction and restitution coefficient. The size of the components is of 
order 10 _1 whereas the magnitude of the next order coefficients in the Fourier Series is much 
smaller. This allows us to truncate the series in Equations (1261) and ([27]) at first order and 
still get good agreement to the actual data, as in Figure HI 

Now that we have a functional form for C(8) and Fb c (6), we can solve for Ti[&) to first 
order in the Fourier components {a c , a' c , aj, a'^}. This gives 

i 

T t (Q) = (F hc )(z-lY^ + ¥(a f sm2e + a' f cos2e) + J2^ i ^ e ~ i ( a c^2e + a' c cos26^j, (28) 

i=i 

where $ and \1/ are variables that depend on the geometry of the force networks and are 
expressible as 

/2tt/3 
dxP c (x) cos(x){l,cos(2x)}. (29) 
-2tt/3 

We can also solve for the constraint in Equation ([25]) . To lowest order in the Fourier 
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FIG. 5: Measured values of the Fourier components from Equation ([2 



components, the constraint equation gives 

T{9) 



$(z - 1) < 1. 



This, combined with Equation (128]) . provides a closed formula for T{. 

J : t {B) = {F hc ) I mm{$(z-l),lY + ¥(z-lY(a f sm29 + a' cos 20) 



(30) 



(31) 



+ ^^(z-iyminf^z-l),!]^ a c sin 29 + a' c cos 26) J, (32) 



to first order in {a c , a' c , a,f, a'y}. 

This solution can now be used to arrive at a constitutive relation for the stress tensor. 
The static stress tensor is given by Equation ([2]), which can be rewritten in two- dimensions 

as 

If I cos 2 6 cos 6 sin 6 

Kp = 77 / deC{9)a{B)F{9) x 

v J \ cos 9 sin 9 sin 9 

where a(9) is the average value of the distance between grains at contact for a given angle. 
In our simulations we observe that a(9) has very little dependence on 9 (of order less than 
10~ 4 ). Thus a{9) = (a). We can also use this same integral form to determine the collisional 
stress tensor by replacing F{9) with F hc {9). 



(33) 
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In this paper we have concentrated on the pressure and shear stress. The pressure is 
given by one-half the trace of Equation (1331) and the shear stress by either off-diagonal 
element, but these two quantities do not fully describe the stress tensor. There is a third 
independent term and, without loss of generality, we use En. Inserting the solution for F{9) 
from Equations (12"UI) and (1521) into Equation (I3"3"j) . we arrive at the following constitutive 
relations that fully describe the stress tensor: 



p bc 



£/€el"l 

min[$(^-l),lf, (34) 
- [a f ¥{z-\Y + a c Y^^\z-l) i TmYi[^(z-l),l} e - i Y (35) 

i=0 J 



p bc 2 ^ 

^P^= £ min[$(z-l),lf + ^_l)^ 

+ y^^-lfmin^-llir L (36) 

i=0 / 

These equations relate the static stress tensor to its collisional values, and properties of 
the force networks. The left hand side (lhs) of each equation gives the difference between 
the static and collisional values of stress. These differences are equal to network properties 
on the right hand side (rhs) of each equation, which are summed over all possible force 
chain path lengths i in the network. A large number of network parameters appear in 
these equations and makes the form rather complicated. This is to be expected, since the 
structure of the force networks is complex and the predictions of Equations (I34H36P span 
system behavior from very dilute granular materials with cf) w to dense granular materials 
with arbitrarily close to <p c (and for all restitution coefficients e). The properties of the 
networks change considerably over this range of parameter space, and they must be retained 
in the Equations. 

Simple scaling relations can be obtained near certain packing fractions. For example, near 
the network transition at tv, Equations (I34H36P predict that (£ig ~ ^ap) ^ ( z — !)• This is 
because, when £/£ e i ~ 1 near u^ c , the deviation from the collisional stress is dominated by 
forces transmitted between nearest neighbors. Therefore the excess number of contacts serves 
as the dominant scaling variable. In contrast, near the jamming transition at u c , networks 
are saturated so that Equation (fT9|) takes its maximum value and all of the kinetic energy 
from each contact is transferred to the network. In this case the stress tensor should depend 
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packing fraction packing fraction 

FIG. 6: (a) The coordination number z and (b) <1> and (Equation (I29|) ) for a wide variety of 
packing fractions and restitution coefficients. The line labels (e-values) in (b) are the same as 
in (a). Although z and depend sensitively on the value of v and e, the value of is always 
approximately 0.83. 

on the size of the networks, and the constitutive relations indeed predict that £lg oc £. This 
scaling is especially interesting since it suggests that the size of the networks is the important 
scaling variable, which might also control other features of the jamming transition. Indeed, 
once £ becomes large, it is clusters of grains, and not individual grains, that serve as the 
basic thermodynamic degrees of freedom. 

Finally, we consider the limit of £/£ e i — > 1 where force networks consist of only two grains. 
In this limit the force network model predicts that E^g = Sj^, which is the appropriate result 
in the dilute regime where kinetic theory holds. 

C. Testing the predictions 

Equations (I34H36P make predictions for all independent components of the stress tensor 
over the complete range of e and for < <f> c . The difference between the static and collisional 
values of stress is related to features of the force networks. These include the anisotropies 
in the contact probability and collisional force {a c , a' c , a/, a^}, the size of the force networks 
£/£ei; the average coordination number z, and a pair of geometric variables {$, \P} that are 
defined in Equation fl29l) and are related to the distribution of contacts on single grains. All 
of these variables except z, $, and \I/ have been measured previously in this paper or in 
Ref. [2|. Our next step is to measure z, and \F 
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In Figure [6] we plot the values of z, $, and ^ as measured in our simulations. We measure 
z by averaging over long-lived contacts, which are pairs of contacting grains that were also 
contacting in the previous time step. This ensures that only the static backbone of the force 
network is considered and that transient contacts do not artificially increase the coordination 
number. This measurement does not depend on the time step, as shown in the inset. We 
also measure $ and as prescribed in Equations (1221) and (|29|) . by averaging over the same 
set of contacts. We observe that \l/ < $ for all granular materials we have considered. 

We have now measured every variable in the constitutive relations of Equations (I34H36I) . 
We can therefore test the validity of the predictions without using any fitting parameters. 
Due to the complexity of the equations, it is convenient to plot the right hand side (rhs) 
of each equation versus the left hand side (lhs). This is shown in Figure [7] using all of the 
data we have collected. Plotted in this way, the data for each component of the stress tensor 
collapses onto the line predicted by the force network model over more than four decades. 
This collapse is especially striking since the variables in the predictions have a wide variance 
as a function of both restitution coefficient and packing fraction. 

The collapse of our data onto the predicted curves suggests that the force network model 
captures an essential property of granular materials over a broad range of densities and 
restitution coefficients. The success of the model is based on visualizing granular materials 
as conglomerates of interacting networks, instead of collections of grains. Thus the packing 
fraction and restitution coefficient, which are grain properties, are substituted by the size, 
coordination, and other properties of the networks. This allows for constitutive relations to 
be determined analytically. 

Finally, it is important to remark that the constitutive relations from the force network 
model have been derived in the limit of perfectly rigid grains, which may not always apply to 
realistic flows with finite grain stiffness. In the case of finite grain stiffness, there is a finite 
speed v c at which forces propagate through the network. Combined with the lifetime of 
the networks r c , this sets a maximum correlation length v c t c , since information can only be 
transferred between a pair of grains if the network exists long enough to propagate it. This 
maximum correlation length is a monotonically increasing function of the grain stiffness. If 
v c t c > £, then the stress tensor can be described by Equations (13411361) . However, if v c t c < £, 
it is necessary to replace the length scale £ with v c r c . Because £ diverges as the material 
approaches the jamming limit and v c t c is always finite, we expect that for a given grain 



22 




10" 3 10~ 2 10" 1 10° 10 1 

rhs 




S -9 -1 (1 1 ? 

10 10 10 10 10 10 

rhs 




10" 3 10" 2 10" 1 10° 10 1 

rhs 



FIG. 7: Tests of the constitutive relations from the force network model, (a) Test of Equation (|34p . 
The left hand side (lhs) of the equation, p ~£ c , is plotted as a function of the right hand side (rhs). 
This collapses the data to the line predicted by the model; (b) The left hand side of Equation ([35]) . 
s ~ b s c , plotted as a function of the right hand side, once again collapsing to the prediction; (c) 
The left hand side of Equation ([36]), 1 ^ bc 11 , plotted as a function of the right hand side. All of 
the plots have been constructed using simulation data for each variable and no fitting parameters 
have been utilized. 
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stiffness there is a critical packing above which v c t c < £. This critical packing fraction is 
always in the inertial regime and must be strictly less than v c . Therefore, very close to 
jamming, the elasticity of grains begins to play an important role. For all other packing 
fractions, the assumption of perfectly rigid grains is valid for natural and experimental flows. 

IV. CONCLUSIONS 

The underlying microscopic interactions between grains have a large influence on macro- 
scopic characteristics of granular flow. We have investigated two models of the stress tensor- 
kinetic theory and the force network model- which make different assumptions about the 
microscopic interactions. Kinetic theory assumes that only binary collisions are relevant 
and calculates the stress tensor based on grain properties, whereas the force network model 
allows simultaneous interaction between many grains and calculates the stress tensor based 
on properties of the resultant force networks. 

For dilute flows, which occur when the size of the force networks is small, kinetic theory 
makes accurate predictions. This is not surprising since small force networks imply localized 
interactions and binary collisions. For dense flows, force networks extend beyond pairs of 
grains and the predictions of kinetic theory no longer match data from simulations. This is 
because grain-grain correlations are induced via the force networks and kinetic theory does 
not take them into account. However, correlations never exist between isolated networks, 
and by constructing the force network model based on network properties, we are able to 
accurately predict the stress tensor for both dilute and dense flows. 

The force network model predicts all independent components of the stress tensor over 
the entire inertial regime and matches data from simulations for more than four decades. 
Extensions of the model could be used to predict other quantities, including the contact 
force distribution function P(f). An integral part of any such theory is the finite size of the 
force networks, which has important effects on the qualitative features of P(f) Further 
extensions could also specify relations between network parameters and thereby simplify 
the constitutive equations (I34H36I) . While it is possible that the relations between network 
parameters are complicated and depend on many factors, it is likely that simple scalings 
exist near the network transition at p^ c . For v w the deviation between the static and 
collisional stress tensors is proportional to (z—1) and it is likely that network parameters also 
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scale with powers of [z — 1). In particular, it would be interesting to probe the dependence 
of £ on (z — 1), although fluctuations in the parameters have complicated our measurement 
of this relation. 

Finally, while we have concentrated on the inertial regime, it is also important to un- 
derstand how natural flows make the transition from dynamics dominated by inertia, to 
quasi-static dynamics, and ultimately to how the system jams. Along this sequence, the 
stiffness of the grains plays an increasingly active role in the dynamics. The force network 
model can accommodate the development of stiffness by incorporating a maximum length 
scale through which forces can propagate. Including this mechanism for a granular material 
moving through v c may help connect the dense inertial regime with the quasi-static regime 
and facilitate a more complete understanding of dry granular materials. 

This work was supported by the William M. Keck Foundation, the MRSEC program 
of NSF under Award No. DMR00-80034, the James S. McDonnell Foundation, the David 
and Lucile Packard Foundation, and NSF Grant Nos. DMR-9813752, PHY99-07949 and 
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